#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
using namespace std;

//Rates
const double Km = 10.0;
const double Kp = 10.0;
const double Gm = 1.0;
const double Gp = 1.0;

//Initial vals
long double G = 1;
long double M = 0;
long double P = 0;

//Time
double timeElapsed = 0.0;

void monteCarlo ();

void monteCarlo() {

        double URN[4];
	
	for(int i=0;i<4;i++){URN[i] = (double)rand()/(double)RAND_MAX;}

	//Propensities
	double KmG = Km * G;
	double GmM = Gm * M;
	double KpM = Kp * M;
	double GpP = Gp * P;

	//Putative Time
	double timeKm;
	double timeGm;
	double timeKp;
	double timeGp;

	//Divide by 0 check
	if (KmG == 0){timeKm = 999999;}else{timeKm = -log(URN[0])/KmG;}
	if (GmM == 0){timeGm = 999999;}else{timeGm = -log(URN[1])/GmM;}
	if (KpM == 0){timeKp = 999999;}else{timeKp = -log(URN[2])/KpM;}
	if (GpP == 0){timeGp = 999999;}else{timeGp = -log(URN[3])/GpP;}

	//Race
	if (timeKm < timeGm && timeKm < timeKp && timeKm < timeGp){
		//Km Wins
		M++;
 		timeElapsed+=timeKm;
        }

	if (timeGm < timeKm && timeGm < timeKp && timeGm < timeGp){
		//Gm Wins
		M--;
 		timeElapsed+=timeGm;
	}

	if (timeKp < timeKm && timeKp < timeGm && timeKp < timeGp){
		//Kp Wins
		P++;
 		timeElapsed+=timeKp;
	}

	if (timeGp < timeKm && timeGp < timeGm && timeGp < timeKp){
		//Gp Wins
		P--;
 		timeElapsed+=timeGp;
	}

}

void print() {
	
	cout << "T = " << setw(6) <<  setprecision (4) << timeElapsed;
	cout << " G = " << setw(3) <<  setprecision (2) << G << " M = ";
	cout << setw(11) <<  setprecision (9) << M << " P = ";
	cout << setw(11) <<  setprecision (9) << P << endl;

}

int main(int argc, char *argv[]){

	srand ( time(NULL) );

	cout << fixed << setiosflags (ios::right);
        
        while(timeElapsed < 300) {
	    print();
	    monteCarlo();
	}
}
